New model for surface fracture induced by dynamical stress 



J0rgen Vitting Andersen* and Laurent J. Lewis^ 
* Department of Mathematics, Imperial College, Huxley Building, 180 Queen's Gate, London SW7 2BZ, England 
* Departement de physique et Groupe de recherche en physique et technologie des couches minces (GCM), Universite de 
Montreal, Case postale 6128, Succursale Centre-Ville, Montreal, Quebec, Canada H3C 3J7 

We introduce a model where an isotropic, dynamically-imposed stress induces fracture in a thin film. 
Using molecular dynamics simulations, we study how the integrated fragment distribution function 
depends on the rate of change and magnitude of the imposed stress, as well as on temperature. 
A mean-field argument shows that the system becomes unstable for a critical value of the stress. 
We find a striking invariance of the distribution of fragments for fixed ratio of temperature and 
rate of change of the stress; the interval over which this invariance holds is determined by the force 
fluctuations at the critical value of the stress. 

PACS numbers: 46.30.Nz, 62.20Mk 



By experience most — if not all — materials will sooner 
or later develop cracks. Yet, a profound understanding 
is largely missing of phenomena such as how cracks initi- 
ate, the formation of networks of cracks and the resulting 
distribution of fragments, the dynamics of crack propa- 
gation, and the collective behavior of many interacting 
cracks. In this Letter we propose a new model that ad- 
dresses, at least in part, some of these questions. In the 
model, an isotropic, dynamically-imposed stress, caused 
by material properties changing in time, induces frac- 
ture in a surface material. The problem is solved using 
molecular-dynamics simulations for a set of beads inter- 
acting with one another via a continuous potential. This 
model should be relevant to many phenomena that are 
known to lead to macroscopic fracture, such as desicca- 
tion jl] Ej or expansion |5|,^| , changes in chemical compo- 
sition [|7|] , changes in temperature J8| , or change of phase 
of the surface layer. 

On the basis of a mean- field argument, we demonstrate 
that the system becomes unstable for a critical value of 
the stress. We find a striking invariance of the distri- 
bution of fragments for a fixed ratio of temperature and 
rate of change of the stress; the interval over which this 
invariance holds is determined by the force fluctuations 
at the critical value of the stress. 

Model — We represent the thin film on a coarse-grained 
scale by beads that mutually interact via a continuous 
potential which we take to be of the Lennard-Jones form, 
4e[(a / r) 12 — (c/r) 6 ], where r = \f\ is the distance between 
two particles. An isotropic stress is imposed by having 
a change in time (t), reflecting a change in the range 
of the interactions on the surface . For simplicity we 
limit ourselves to the case where the material is initially 
unstressed and o~(t) decreases monotonically with time. 
This corresponds to a surface where the induced stress 
makes the material rupture in a state of tension. 

The dynamics of the beads obeys Newton's second 
equation, i.e., the system is simulated using molecular 



dynamics (MD). We assume the surface layer to be in 
contact with a heat bath at temperature T; this is done 
by periodically rescaling the velocities to a fixed kinetic 
energy p0[ . The units are chosen so that the mass 
m = e = 1. In its initial (stress-free) state, the surface 
layer consists of a triangular lattice with lattice constant 
ao = 26<7o, where oq = a (to). Periodic boundary con- 
ditions are used to eliminate surface effects. The conse- 
quence of decreasing o~(t) is to put all beads under tensile 
stress, i.e., each bead feels attracted by its neighbors. We 
assume a(t) to decrease linearly in time until it attains a 
final value cry at time tf, whereafter it remains constant. 
An effective strain parameter of the overlayer is defined 
by s(t) = [do — a(t)]/o- . The rate of change ("speed") 
of a is denoted v = da(t)/dt. 

We are interested in the fracture pattern at t = oo 
which is obtained in practice by choosing a large enough 
tf, whose value depends on s(tf), v, and T; the latter 
three parameters determine completely the fracture pat- 
tern. In order to calculate the probability P(f) for having 
a fragment of size /, we discretize the system into cells 
of size cr(t). A fragment is then defined as a cluster of 
beads that are nearest or next-nearest neighbors to one 
another. 

As a changes with time, each bead will evolve from a 
position of global energy minimum to a local minimum 
state. The local minimum energy state is stable, how- 
ever, for cr(t) close to 00, since the system would need 
instant cooperative motion of all the beads in order to 
rearrange into the global minimum-energy state whose 
lattice parameter is a = 2ea(t). Due to the many body 
nature of the system, each bead will see an energy land- 
scape that changes as the positions of neighboring beads 
change, and as a changes in time. The cooperative mo- 
tion of the beads create dynamical and spatial barriers 
between, on the one hand, local metastablc minimum en- 
ergy states, and, on the other hand, the global minimum 
energy state. 
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For increasing values of a, the initial configuration 
eventually become unstable. Neglecting fluctuations in 
the positions of the beads, each will experience a mean 
field potential from its nearest neighbors given by: 

V(r, a) = l2e[C-r - C~f + (^) 12 - i^f]- 
r r la — r la — r 

We are interested in the behavior of V(r, a) at the point 
r = a + 5 with S small. Expanding the above to fourth 
order in <5, we find: 

V(S, a) = 12e(-) 6 {2[(-) 6 - 1] + [156(-) 6 - 42](-) 2 
a a a a 

+ [32760(-) 6 - 3024](-) 4 + 0((-) 6 )}. 
a a a 

Thus, for 5 small, the potential seen by a bead changes 
from a harmonic single- well to a double-well potential as 
a decreases. This happens when V {S,&)\s=o changes 
sign, that is for a c = (7/26) 1 / 6 a = (7/13) 1 /6 a - « 
0.90fio- In general, the existence of a critical a c for an ar- 
bitrary interaction V(r, a) is equivalent to V (r, a)\ r=a = 
having a solution. As a(t) approaches <r c from below, 
one large fluctuation eventually takes place bringing one 
of the beads close to it's new local minimum energy po- 
sition. A cascade of similar events then spreads out from 
beads adjacent to that which first broke the configura- 
tional symmetry. The extent of the propagation of this 
cascade of events, and the subsequent fracturing of the 
system, depends, as we will see, on s(tf), T, and v, as 
well as on the fluctuations of forces when a = a c . 

Results — Figs, la-c show snapshots of one system 
for different values of stress but fixed temperature and 
stress speed. Fig. la corresponds to a stress a(t) slightly 
larger than a c . The very first cracks have appeared and 
shortly after the system completely disintegrates into 
many pieces, characterized by a macroscopic Young's 
modulus that goes to pp|]. This has happened in 
Fig. lb. In Fig. lc, we have the final state of the system 
when the stress no longer varies in time. The effect of 
varying the speed v can be seen in Figs. Id and le: here, 
the initial conditions are the same as in Figs, la-c, but v 
is 8 times smaller. The stress in Fig. Id is the same as 
in Fig. lb; clearly, a smaller rate of change of the stress 
gives the system longer time to respond so that the posi- 
tions of the beads are correlated over a longer distances 
and the cracks are straighter. As a result, the fragments 
in the final configuration, Fig. le, are larger than they 
are under a rapidly- varying stress (compare Fig. lc). 

If T = the absence of thermal fluctuations would 
mean that the system remains in its initial state and 
never breaks, despite the fact that the energy difference 
between initial and stressed states increases as a(t) de- 
creases. For T ^ |l2j and v — > oo, on the other hand, 
the rupture of the system is completely dominated by 
fluctuations, in which case the probability density P(f) 



for having a fragment of given size / is given by a bino- 
mial distribution P(f) = Kr§ ,/)(^K(§) 6_ ^\ since each of 
the 6 neighbors of a given bead has probability | of form- 
ing a cluster with that bead. For finite (T, v), finally, the 
fracturing is determined by the coherent motion of the N 
beads. In Fig. 2a-b we show the cumulative probability 
distribution P>(/) for a given T and different v; as we 
have seen above, the smaller the value of v, the larger the 
fragments. In Fig. 2a, s(tf) = 0.5, whereas s(tf) = 0.75 
in Fig. 2b. In order to calculate P>(f), we have averaged 
over 200-500 N = 100 systems with different initial con- 
figurations, all at the same temperature T. (We chose 
to use many small systems rather than few large ones 
in order to get better statistics). Finite-size scaling of 
P>(f) is shown in the inset of Fig. 2a, which allows us 
to extend our results (for the given (T, v)) to the case 
N — ► oo. The lines are fits to a log- normal distribution; 
clearly, the data suggest this form of P>(/) for large v. 
This is the signature of a fracturing process that happens 
in a multiplicative manner |T^| , where a given piece at a 
random point breaks into two pieces, which themselves 
randomly break into two other pieces, etc. For very small 
v, P>(f) crosses over to a Heaviside theta function, since 
in this case breakdown happens due to one large crack 
spanning the whole system. The speed for which P>(/) 
can no longer be described by a log-normal distribution 
depends on T and N, and is due to finite size effects. 

An instantaneous change in a means a change in both 
the magnitude and the fluctuations of the forces. We find 
the system to respond in a qualitatively different manner 
to changes in a depending if it is < a c or > a c . For a 
broad range of speeds v, we find the average magnitude 
of the force on the beads, F = \fi\/N, and its fluc- 
tuation, SF = Y^,i y/ff — F 2 /N, to be independent of v 
for a(t) < a c . We have also calculated the characteristic 
length, of the stress field F[r(t)] by taking the first 
moment of the radial averaged structure factor S(k, t). In 
0, a coarsening phenomenon of F[r(t)] prior to the first 
fracture was found to be crucial for the subsequent rup- 
ture of the system; in the present model, we observe no 
time evolution of for <r(t) < a c , and £(t) ~ a. (how- 
ever, when the first macro cracks appear, £(i) increases 
dramatically). Therefore, the observed dependence of 
P>(f) on v must be due to the way the system responds 
to changes in a after the a c point has been passed. 

Whether or not the system has time to counteract 
the imposed stress passed a c depends on the timescale 
over which changes in a take place compared to the re- 
sponse time of the system; the latter is determined by 
the random thermal motion, i.e., kinetic energy Ek, of 
the beads. The ratio of these two timescales is thus given 
by k = av^ 1 / (mi a / E^ 2 ) — y/E^/v. One therefore ex- 
pect systems with the same value of n to fracture in the 
same way. The fracture is expected to be dominated by 
fluctuations for k <C 1, whereas for k ^S> 1 it will have 
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time to respond to the changing stress in a correlated 
manner. This is in fact verified in Fig. 3 which shows a 
remarkable invariance of P>(f) over almost 3 decades in 
temperature for systems with two different values of k. 
The lowest and highest temperature in Fig. 3 for which 
the invariance of P>(f) no longer holds, and the subtle 
temperature dependence at intermediate values, can be 
understood from the dependence on stress of the force 
fluctuations SF, shown in Fig. 4 for the same values of T 
as Fig. 3. Because of fluctuations, different temperatures 
lead to a critical s c (defined as the s for which 8F has 
its minimum) slightly different from the mean field value 
of s c = (uo — cr c )/(jQ = 0.10. The small temperature 
dependence of P>(f) at intermediate temperatures can 
then be understood in terms of a slight increase of SF(s c ) 
with T, since one would expect larger force fluctuations 
at s c to lead to smaller fragments. As seen in Fig. 4, 
the only exception to this is the case of the highest T 
where, on the contrary, a large SF(s c ) leads to a large- 
fragment tail in P> (/). The reason for this is that T is so 
high that coalescence of already-formed fragments takes 
place; coalescence is not observed for lower T. Finally 
one also notes from Fig. 3 that deviations in P>(f) oc- 
cur for very low temperatures, where the simple scaling 
argument leading to invariance of P>{f) under a given k 
apparently no longer holds. 

Conclusion - We have introduced a model where 
a dynamically-imposed stress induces fracture in a thin 
film. Using molecular-dynamics simulations, we have 
shown the accumulated fragment distribution function to 
obey a log-normal distribution characteristic of fractur- 
ing processes which happen in a random multiplicative 
manner. A mean field argument shows how the system 
undergoes an instability for a critical value of the imposed 
stress. We find a striking invariance of the fragment dis- 
tribution function for a given ratio of temperature and 
speed of stress; the interval over which this invariance 
holds, is determined by the force fluctuations at the crit- 
ical value of the stress. 
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FIG. 1. Snapshots of a N = 1600 system at different times t, with different change of strain rate v and different final 
strain s(£/). The initial configuration is the same, and T = 6.25 x 10~ 5 , in all cases, (a) s(t) — 0.14, v = 0.0125, (b) 
s(t) = 0.25, v = 0.0125, (c) s(t = t f ) = 0.5, v = 0.0125, (d) s(t) = 0.25, v = 0.0015625 and (e) s(t = t f ) = 0.5, v = 0.0015625. 



FIG. 2. Cumulative probability distribution P>(f) for finding a given fragment of area larger than /, and T = 6.25 x 10 5 
(a) s(t f ) = 0.5 ; v = 0.025 (o), 0.0125 (+), 0.00625 (□), 0.0042 (x), 0.003125 (A) and 0.0015625 (*). The lines to are fits to a 
log-normal distribution. Inset: finite size scaling with s(t f ) = 0.5, v = 0.0125; N = 100 (x), 400 (A) and 900 (*); s(t f ) = 0.75, 
v = 0.0125; N = 100 (o), 400 (+) and 900 (□). (b) s(t/) = 0.75 ; v = 0.0375 (o), 0.01875 (+), 0.009375 (□), 0.0046875 (x), 
0.003125 (A), 0.002679 (*) and 0.002344 (small white circle). The lines to are fits to a log-normal distribution. 



FIG. 2. Inset to Fig. 2a. 



FIG. 2. Fig.2b 



FIG. 3. P>(/) versus / for fixed value of k = E\J 2 /v and s(t f ) for a N = 100 system, (i) k = 1.03, s(t f ) = 0.75, and 
(T,v) = (6.4 x 10" 2 , 0.30) (o), (1.6 x 10~ 2 , 0.15) (+), (4 x 10~ 3 , 0.075) (□), (10~ 3 , 0.0375) (x), (2.5 xl0~ 4 , 0.01875) (A) and 
(6.25 xl0~ 5 , 0.009375) (*). (ii) k = 1.55, s(t f ) = 0.5, and (T,t>) = (6.4 x 10~ 2 , 0.20) (large black circle), (1.6 x 10~ 2 , 0.10) 
(black circle), (4 x 10~ 3 , 0.05) (small white circle), (10~ 3 , 0.025) (white circle), (2.5 xl0~ 4 , 0.0125) (large white circle) and 
(6.25 xlO -5 , 0.00625) (small black circle). 



FIG. 4. SF versus s for n = 1.55 and (T,v) = (6.4 x 10" 3 , 0.20) (o), (1.6 x 10" 2 , 0.10) (+), (4 x 10~ 3 , 0.05) (□), (10~ 3 , 
0.025) (x), (2.5 xl0~ 4 , 0.0125) (A) and (6.25 xlO -5 , 0.00625) (*). 



FIG. 4. Fig.4-inset 
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